Testing quantum electrodynamics in the lowest singlet states of beryllium atom 
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High-precision results are reported for the energy levels of 2 1 5' and 2 1 P states of the beryllium atom. Calcu- 
lations are performed using fully correlated Gaussian basis sets and taking into account the relativistic, quantum 
electrodynamics (QED), and finite nuclear mass effects. Theoretical predictions for the ionization potential of 
the beryllium ground state 75 192.699(7) cm -1 and the 2 1 P ->• 2 1 S transition energy 42 565.441(11) cm" 1 
are compared to the known but less accurate experimental values. The accuracy of the four-electron computa- 
tions approaches that achieved for the three-electron atoms, which enables determination of the nuclear charge 
radii and precision tests of QED. 

PACS numbers: 31.30.J-, 31.15.ac, 32.10.Hq 



Spectroscopic standards for the energy transitions of the 
beryllium atom have been established many years ago (1962) 
in the experiments by Johansson. In most cases an accuracy 
of 0.01 -0.02 A fl] has been reached. The only more precise 
beryllium energy level determination comes from the exper- 
iment by Bozman et al. [2] in which the transition energy 
of 21 978.925 cm" 1 between the ground and the 2 3 Pi state 
accurate to about 0.01 cm -1 (0.002 A) has been measured. 
Seaton, through a fit to a collection of excited states data, has 
determined the ionization potential (IP) of the ground state to 
be 75 192.56(10) cm -1 y[]. Later, this quantity has been im- 
proved by Beigang et al., who obtained 75 192.64(6) cm" 1 
0]. Contemporary high precision calculations ^ are in 
good agreement with the rather old experimental IP values. 
Nevertheless, the precision of the data available for beryllium 
is far from being satisfactory compared to the exquisite accu- 
racy of the modern atomic spectroscopy. The level of the ab- 
solute precision achieved in modern measurements for three - 
electron systems is as many as four orders of magnitude 
higher than that obtained in the case of beryllium. As it has 
been shown for two- and three-electron atoms, the availability 
of such accurate data, in connection with good understand- 
ing of the underlying atomic theory, opens up access to such 
interesting applications like the determination of the nuclear 
charge radius or precision tests of the quantum electrodynam- 
ics (QED). 

Remarkable advances in theoretical methods make it pos- 
sible to approach the spectroscopic accuracy for the energies 
and transition frequencies of few-electron atoms. This chal- 
lenge requires precise treatment of the electron correlations as 
well as inclusion of relativistic and quantum electrodynamic 
effects. The concise approach, which accounts for all the ef- 
fects beyond the nonrelativistic approximation, is based on the 
expansion of the energy levels in the fine structure constant a 
(see Eq. (Q~|i below). This method has been successfully ap- 
plied in recent years to light atomic and molecular systems 
Il9l 4l2n . The frontiers in this field of research have been estab- 
lished by the calculation of higher order (m a 6 ) corrections to 



helium energy levels [13] and of ma 7 corrections to helium 
fine structure ifTill . 

Up to now, the precision of theoretical predictions for the 
beryllium states with the non-vanishing angular momentum 
has been severely limited by the accuracy of the lowest-order 
relativistic (m a 4 ) and QED (m a 5 ) corrections. The most ac- 
curate calculations including the relativistic corrections was 
performed 20 years ago by Chung and Zhu [15] at the full- 
core plus correlation level of theory, whereas the QED effects 
have merely been approximated from hydrogenic formulas 
lfl6i[l7tl . This approach has turned out to be unsatisfactory for 
it has led to a significant discrepancy between theoretical pre- 
dictions and experimental excitation energies. For instance, 
the theoretical result IU5I1 is by 3.45 cm" 1 higher than the 
experimental value 2 349.329(10) A of the 2 1 P -> 2 1 S' tran- 
sition energy measured with 0.18 cm -1 uncertainty |Q]]. The 
fact that the theoretical excitation energy is higher than the 
experimental value may indicate that correlation effects have 
not been incorporated satisfactorily. Such a disagreement can 
only be resolved in an unequivocally more accurate computa- 
tion of nonrelativistic energies as well as the relativistic and 
QED effects using the explicitly correlated wave functions. 
Recent nonrelativistic calculations of low-lying P- and In- 
states with the relative precision of an order of 10" 10 -10 
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[18 



1911 represent a step in this direction. 
In this paper, we present the first complete and highly ac- 
curate treatment of the leading relativistic (ma 4 ) and QED 
(m a 5 ) effects for a four electron atomic P-state. We also im- 
proved results for the 2 1 S f state, which permitted us to push 
the accuracy of the theoretical predictions of the 2 X P — > 2 1 S 
transition energy beyond the experimental uncertainty. Addi- 
tionally, in combination with the previously reported very ac- 
curate data on beryllium cation 111 Oil , we obtained an improved 
ionization potential with the accuracy an order of magnitude 
higher than that of the available experimental values. 

In our approach, we expand the total energy not only in the 
fine structure constant a s» 1/137 but also in the ratio of the 
reduced electron mass to the nuclear mass ?y = —/.i/tun — 
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— m/(m + mjv) » 1/16424. This way we reduce the isotope 
dependence to the prefactors only. In terms of these two pa- 
rameters, the energy levels can be represented as the following 
expansion 

E= ma 2 + V 8^] + ma 4 £^ 

£(5,0) + 

ma £<■*& + .... (1) 

Each dimensionless coefficient £(™> n ) is calculated separately 
as an expectation value of the corresponding operator with 
the nonrelativistic wave function. The leading contribution 
£( 2 .o) = g Q j s an eigenvalue of the Schrodinger equation with 
the clamped nucleus Hamiltonian 



Ho* = So*, Ho 



The key to obtaining high-precision results is the use of a very 
accurate trial wave function '5, which contains all the inter 
particle distances explicitly incorporated. We express f as a 
linear combination of N four-electron basis functions ipi 



N 



tpi = A[(t>i(ri, r 2 , r 3 , r 4 ) x] , (3) 



where A is the antisymmetry projector, \ = 
\ (ti4-a - 4-ita) (ts-U - 4-3T4) is the singlet spin func- 
tion constructed using electron spinors. The spatial function 
<f> is the explicitly correlated Gaussian (ECG) function for S- 
and P-state, respectively 

4> s = exp [ - Y w a r\ - ^ u ab rl b ] , (4) 

a a<b 

$ P = r x exp [ - Y w a r\ - ^ u ab r 2 ab ] . (5) 



a<b 



The main advantage of these Gaussian functions is the avail- 
ability of analytical forms of the integrals required for matrix 
elements of the Hamiltonian Ho 



f(ni, ■ 
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a<b 



Among all the integrals represented by the above formula we 
can distinguish two subsets used in our calculations. The first 
subset contains the "regular" integrals with the non-negative 
even integers rij such that m < fli, where the shell pa- 
rameter CI i = 0, 2, 4, ... . The second subset permits a sin- 
gle odd index n.i > —1 for which J2i n i — ^l/r i^i/r — 
— 1,1,3,...) and is related to the components of the Coulomb 
potential. To systematize the use of the ECG basis sets we 
re-derived the recurrence scheme for the generation of both 
the classes of integrals from the master expression 1 20l 21 1 . 
An advantage of such approach is the possibility of a gradual 



extension of calculation to the states with higher angular mo- 
menta. The sets of integrals employed in a specific case can 
be characterized using the fl shell parameters. For instance, 
the matrix elements of the nonrelativistic Hamiltonian require 
integrals with fij = 2, fii/ r = —1 for 5-states (Eq. (0]i), and 
fli = 4, Q,x/ r = 1 for P-states (Eq. © ). If, additionally, gra- 
dients with respect to the nonlinear parameters are to be used, 
both shell parameters have to be increased by two. 

To control the uncertainty of our results we performed 
the calculations with several basis sets successively increas- 
ing their size by a factor of two. From the analysis of 
convergence we obtained the extrapolated nonrelativistic en- 
ergies and mean values of the operators presented in Ta- 
ble U The largest wave functions optimized variationally 
were composed of 4096 and 6144 terms for the S- and 
P-state, respectively, leading to the nonrelativistic energies 
gi 2 '°)(2 1 S) = -14.667356494 9 a.u. and £( 2 -°)(2 1 P) = 
— 14.473 451 33 4 a.u. These upper bounds improve slightly 
those obtained by Adamowicz et al. 

The other coefficients of expansion (fl]i are calculated as 
mean values with the nonrelativistic wave function The 
nonrelativistic finite mass correction is given by £^ 2 ^ = 
£( 2 >°) _ ^ a<fe {p a ■ Pb}- In order to calculate the leading rel- 
ativistic corrections £ ' 4, °) = {H 1 -- 4 ' ^) we consider the Breit- 
Pauli Hamiltonian lE2tl . which for the states with vanishing 
spin can be effectively replaced by the form 



H (4,0) = 



7r Z a 



a<b 



Tr$ 3 (r ab 



;P a 



(7) 



Tab 
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' ab 



Since the ECG basis does not reproduce the cusps of the wave 
function, a slow convergence becomes evident for relativistic 
matrix elements of the Dirac S and the kinetic energy opera- 
tor p\. To speed up the convergence, the singular operators 
can be transformed into their equivalent forms, whose behav- 
ior is less sensitive to the local properties of the wave func- 
tion. For the Dirac 6 expectation value, such a prescription 
has been proposed by Drachman |23|]. For example, from di- 
rect calculation with the basis size of 4096 for S'-state, we 
get (6(r a )) = 35.366 89 . . ., while using the Drachman regu- 
larization approach we improve the convergence by three or- 
ders of magnitude (see Table Regularization methods have 
also been applied for the beryllium ground state in the for- 
mer paper [5], nonetheless the present results are more ac- 
curate by 2 orders of magnitude due to the better optimized 
wave function. For P-states, the expectation values of the 
relativistic and QED operators as well as of the Bethe log- 
arithm have been unavailable in literature to date. Analo- 
gous calculations of relativistic terms in the ECG basis have 
been performed only for P-states of the four-body positron- 
ium molecule M24I1 . Methods for evaluation of additional in- 
tegrals "1/r 2 " and "l/(r a rb)" of the form © required for 
the regularized operators of the Breit-Pauli Hamiltonian have 
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been developed, resulting in computationally tractable recur- 
sive expressions derived from corresponding master integrals. 
These have been presented in the original paper only for three- 
body systems 12111 . 



TABLE I: Expectation values of various operators with nonrelativis- 
tic wave function for 2 1 S and 2 1 P states of beryllium atom (in a.u.). 



Operator 



2 1 S 



2 X P 



Pa\ 



Ho 

Pa ■ Pb 

Pi 

S(r a ) 

S(r ab ) 

j j-llhl 
-H 



1 )pI 



P{r~J) 
In ko 



-14.667 356 498(3) 
0.460 224112(8) 
2165.6301(9) 
35.369 002 6(6) 
1.605 305 33(9) 
1.783 64819(15) 
-7.326 766(3) 
5.750 46(2) 



-14.473 45137(4) 

0.434 81125(13) 
2133.321 1(12) 
34.897 914 6(8) 
1.567 943 6(2) 

1.624185 8(5) 
-7.09715(8) 
5.752 32(8) 



The calculation of the leading QED corrections is the main 
challenge of this work. It is particularly laborious because we 
deal with the states of the non-vanishing angular momentum. 
The explicit form of the m a 5 terms is given by J25 , 26 ] 



t ~ 3 



19 

— +ln(a- 2 ) -lnfc 
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In a 



a 



a<b 



a<b 



This expression contains two highly nontrivial terms: the 
Bethe logarithm In ko and the so-called Araki-Sucher distri- 
bution P{r~i). In ECG basis, the latter exhibits exception- 
ally slow convergence when evaluated directly from its def- 
inition. The regularization is in this case mandatory if one 
aims at a high accuracy of the final results. For this purpose 
we extended the original Drachman's idea and obtained the 
following regularized form for the distribution [E7II 



(Pa— Pc 



Tab 



(9) 



+ {^{\+ 1 )8{r ab ) + 2{£ -V) 



In r ab 

Tab 



As we can see, new classes of the integrals containing factors 
of the form "In r/r", "In r/r 2 ", "lnr a / (r a rb)" arise. With the 
master integral, such integrals can be expressed analytically in 
terms of elementary and Clausen functions. 

The evaluation of the Bethe logarithm is the most time con- 
suming part of the calculations. Formulas for such calcula- 
tions with the ECG functions have been presented in the for- 
mer work devoted to lithium atom [28] and later on applied 
to the beryllium ground state fl. In principle, we use the 
same compact integral representation of Bethe logarithm as 
before, however, for the 2 1 P state such calculations become 
more sophisticated. Compared to the ground state, which 



through the momentum operator is coupled only with the vir- 
tual 1 P states, the 2 X P state requires a complete set of the X S, 
l P e , and 1 D intermediate states. These three types of states 



can be well represented in the bases 4>s, e 



j>3 



4>s, and 



((r* r J b + r J a r£)/2 — 1/3 r\ ri) cf>s, respectively. Evalua- 
tion of f(t)/w in the limit of to = is clearly established nu- 
merically from the Thomas-Reiche-Kuhn sum rule for dipole 
oscillator strengths (P (U - So)- 1 P) = 3 Z/2. This value 
is useful in judging the completeness of the virtual states and 
estimation of uncertainties. 

Because of principal difficulties, the ma 6 corrections in 
their full form were evaluated only for two electron atoms 
lfl3Tl . Therefore, for the four-electron beryllium atom we use 
the following approximate formula based on the hydrogen 
atom theory B29I1 



£(6,0) 



427 
"96" 



21n(2) 



This approximation includes the dominating electron-nucleus 
one-loop radiative correction and neglects the two-loop radia- 
tive, electron-electron radiative, and the higher order relativis- 
tic corrections. On the basis of the experience gained in he- 
lium calculations lfl3ll . we estimate, considering higher charge 
of the beryllium nucleus, that the neglected terms contribute 
less than 20% to the overall m a 6 correction. 



TABLE II: Components of the 2 l P — 2 1 S transition energy and 
the ionization potential (IP) for 9 Be atom in cm -1 . COD ATA l3lll 
inverse fine structure constant a -1 = 137.035 999 074(44) and the 
nuclear mass mjv( 9 Be) = 9.012 182 20(43) u (3(3]. 



Operator 


2 1 P-2 1 S 


IP(2 1 5') 


m a 


42 557.255(6) 


75 190.543(4) 


2 

ma r\ 


-2.930 72 


-4.675 65 


4 

m a 


12.167(1) 


7.414 0(8) 


m a 5 


-1.003 3(14) 


-0.557 7(3) 


m a 6 


-0.045(9) 


-0.025(5) 


Total 


42 565.441(11) 


75 192.699(7) 


Theory Q5J 


42 568.80 




Theory [6] 




75 192.667(19) 


Theory [5] 




75 192.510(80) 


Experiment [1] 


42 565.35(18) 




Experiment [3] 




75 192.50(10) 


Experiment [4] 




75 192.64(6) 



Except for all the coefficients of the expansion (Q]) 

are evaluated in complete, i.e. no approximation is introduced 
nor any physical effect of given order is omitted. Therefore, 
the uncertainties given in Table HI] refer to the incomplete- 
ness of the basis set. On the basis of our former work II 1 1 
on Be + , the higher order in a and -q contributions, namely 



m a r\ and m a rj are estimated to be less than 0.001 cm 
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to both the transition energy and ionization potential, and thus 
they are negligible when compared to the present uncertainty 
of 0.01 cm -1 . We note in passing that in Table ITT1 for the val- 
ues without the uncertainty all the quoted digits are certain. 
In evaluation of the IP value we used the ground state energy 
level of the beryllium cation £(Be+) = -14.325 836 7 a.u. 
calculated with Hylleraas wave functions 11 Oil . 

The accuracy of 0.011 cm -1 for the transition 2 1 P — 2 1 S 
and 0.007 cm -1 for the ionization energy has been achieved 
due to the recent progress made in two directions. The first 
one, essential to reach this accuracy, is the improvement in 
the optimization process of the nonrelativistic wave functions 
leading to the overall increase in numerical precision. The 
second direction is the complete treatment of the leading rel- 
ativistic and QED effects. More specifically, the approach to 
effectively calculate the many electron Bethe logarithm and 
mean values of singular operators, like the Araki-Sucher term, 
has been developed (53. Particularly, an extension of the nu- 
merical methods for relativistic and QED corrections on P- 
states of a four-electron system is presented here. 

The ma 5 and ma 6 terms involve the interaction of the 
electrons with the vacuum fluctuations of electromagnetic 
field, the electron-positron virtual pair creation, and the re- 
tardation of the electron-electron interaction. Results of Ta- 
ble HI] show clearly that taking into account such energeti- 
cally subtle QED effects is unavoidable in to reach the agree- 
ment between the experiment and the theory and that it en- 
ables testing of QED. For example, the overall contribution of 
the QED effects to the 2 1 P — 2 1 S transition energy amounts 
to 1.048(9) cm -1 and is an order of magnitude higher than 
the experimental uncertainty. Currently, the accuracy reached 
by theory for the transition frequencies exceeds by an order 
of magnitude that of the known measurements for beryllium 
atom. At this level of accuracy we are able to resolve the 
2 l P — 2 1 S line discrepancy of 3.45 cm -1 between the ex- 
periment lfl5ll and theory in favor of the latter. Although, the 
available semi-empirical results of the ionization energy [4] 
agree well with the more accurate theoretical results obtained 
here, we hope that the increased level of accuracy of the the- 
oretical predictions will be a stimulus for new, more accurate 
measurements. 

Uncertainty of our results comes mainly from the neglect 
of the higher order relativistic and QED corrections of the or- 
der a 6 ' 7 . Evaluation of these term sets the direction of our 
future efforts. Also the numerical accuracy of the nonrela- 
tivistic energy has to be improved to achieve further progress 
in theoretical predictions. 

The recursive method of evaluation of the integrals (f)J em- 
ployed in this work allows an application of the ECG func- 
tions to the states with non-vanishing angular momentum. 
It establishes a framework for a high accuracy studies of 
the fine structure and the hyperfine splitting in the beryllium 
atom. The isotope mass shifts can also be precisely calcu- 
lated. However, an accurate experimental data are necessary 
to enable an extraction of a nuclear-model-independent charge 
radii from isotope shifts by combining high-accuracy mea- 



surements with atomic theory. This is of special interest for 
the halo nuclei (e.g. n Be) for which analogous results ob- 
tained recently from the 2 2 P — 2 2 S transition in beryllium 
cation require a confirmation. Systematic extension to 

transition energies involving instates is mandatory to resolve 
other severe discrepancies between theory and measurements 
pointed out by Chung and Zhu IU5I1 like e.g. the largest one 
of 7.38 cm -1 for the 3 1 D — 2 X S transition. To our knowl- 
edge, no precise calculations of the Bethe logarithm with such 
nontrivial angular momentum structures have been performed 
for any few-electron systems so far (even for helium). The 
methodology presented in this paper opens up a route towards 
removing such obstacles and, what is more, is very promising 
in applications to five- and six-electron systems. 
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